function [label_lkg, Z, K, g] = LKGr(K, alpha, beta, gamma, mu, nCluster)
% initialization
[nn, nn, m] = size(K);

maxIter=100;
Y1=zeros(nn);
Y2=zeros(nn);
J=eye(nn);
g = ones(1,m)/m;
A = zeros(1,m);
B = ones(1,m);
C = 1;
D = zeros(nn);
lb =zeros(1,m);
% Z=rand(n);
for j = 1:m
    D = D + g(j)*K(:,:,j);
end

M =zeros(m);
for a = 1:m
    for b = 1:m
        M(a,b) = trace( K(:,:,a)*K(:,:,b));
    end
end
W =D;
%main function
for iter = 1:maxIter
    %**********************************************
    % Update Z
    %**********************************************
    if exist('Zold','var')
        Zold = Z;
    end
    Z=(D+mu*eye(nn))\(mu*J+Y1+D);
    Z = Z .* (Z > 0);
    if sum(sum(isnan(Z))) + sum(sum(isinf(Z))) > 0
        L=(Zold+Zold')/2;
        V = spectral_clustering(L, nCluster);
        label_lkg = litekmeans(V, nCluster, 'MaxIter', 100, 'Replicates', 50);
        return;
    end
    %**********************************************
    % Update K
    %**********************************************
    HI=zeros(nn);
    for j = 1:m
        HI = HI + g(j)*K(:,:,j);
    end
    Dold=D;
    %D=(mu*eye(nn)+2*gamma*eye(nn))\(mu*W+Y2-eye(nn)*0.5+Z'-0.5*Z*Z'+2*gamma*HI);
    D = (mu*W+Y2-eye(nn)*0.5+.5*(Z+Z')-0.5*(Z*Z')+2*gamma*HI)/(mu + 2 * gamma);
    D = D .* (D > 0);
    if sum(sum(isnan(D))) + sum(sum(isinf(D))) > 0
        L=(Z+Z')/2;
        V = spectral_clustering(L, nCluster);
        label_lkg = litekmeans(V, nCluster, 'MaxIter', 100, 'Replicates', 50);
        return;
    end
    %**********************************************
    % Update J
    %**********************************************
    G=Z-Y1/mu;
    [U, X, V] = svd(G, 'econ');
    diagX = diag(X);
    svp = length(find(diagX > alpha/mu));
    diagX = max(0,diagX - alpha/mu);
    if svp < 0.5 %svp = 0
        svp = 1;
    end
    J= U(:,1:svp)*diag(diagX(1:svp))*V(:,1:svp)';
    J = J .* (J > 0);
    %**********************************************
    % Update W
    %**********************************************
    H=D-Y2/mu;
    [U, S, V] = svd(H, 'econ');
    diagS = diag(S);
    svp = length(find(diagS > beta/mu));
    diagS = max(0,diagS - beta/mu);
    if svp < 0.5 %svp = 0
        svp = 1;
    end
    W= U(:,1:svp)*diag(diagS(1:svp))*V(:,1:svp)';
    W = W .* (W > 0);
    %**********************************************
    % Update g
    %**********************************************
    for a = 1:m
        A(a)= 2*gamma*trace( D*K(:,:,a));
    end
    options = optimoptions('quadprog','Algorithm','interior-point-convex','Display','none');
    g = quadprog( 2*M*gamma,-A,[],[],B,C,lb,[],[],options);
    
    Y1=Y1+mu*(J - Z);
    Y2=Y2+mu*(W - D);
    
    mu=mu*1.1;
    
    if((iter>5)&&(norm(D-Dold,'fro') < norm(Dold,'fro') * 1e-3))
        break
    end
end

L=(Z+Z')/2;
V = spectral_clustering(L, nCluster);
label_lkg = litekmeans(V, nCluster, 'MaxIter', 100, 'Replicates', 50);
end
function [ V] = spectral_clustering(W, k)

D = diag(1./(eps+sqrt(sum(W, 2))));
W = D * W * D;
[U, s, V] = svd(W);
V = U(:, 1 : k);
V = normr(V);

%ids = kmeans(V, k, 'emptyaction', 'singleton', 'replicates', 100, 'display', 'off');
% ids = kmeans(V, k, 'start','sample','maxiter', 1000,'replicates',100,'EmptyAction','singleton');
end
